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ABSTRACT 



' Context. For magnetically driven events, the magnetic energy of the system is the prime energy 

reservoir that fuels the dynamical evolution. In the solar context, the free energy {i.e., the energy 
in excess of the potential field energy) is one of the main indicators used in space weather fore- 



C/3 ' casts to predict the eruptivity of active regions. A trustworthy estimation of the magnetic energy 

■ 

is therefore needed in three-dimensional (3D) models of the solar atmosphere, e.g., in coronal 
fields reconstructions or numerical simulations. 

■ Aims. The expression of the energy of a system as the sum of its potential energy and its free 

numerical realizations on a discrete grid, this property may be only approximately fulfilled. We 

cn 
o 

CO ' the amount of free energy present in a magnetic configuration. 



energy (Thomson's theorem) is strictly valid when the magnetic field is exactly solenoidal. For 
numerical realizations on a discrete grid, this property may be only approximately fulfilled. We 
show that the imperfect solenoidality induces terms in the energy that can lead to misinterpreting 



Methods. We consider a decomposition of the energy in solenoidal and nonsolenoidal parts 
, which allows the unambiguous estimation of the nonsolenoidal contribution to the energy. We 

■ apply this decomposition to six typical cases broadly used in solar physics. We quantify to what 



extent the Thomson theorem is not satisfied when approximately solenoidal fields are used. 

Results. The quantified errors on energy vary from negligible to significant errors, depending on 
the extent of the nonsolenoidal component of the field. We identify the main source of errors and 
analyze the implications of adding a variable amount of divergence to various solenoidal fields. 
Finally, we present pathological unphysical situations where the estimated free energy would 
appear to be negative, as found in some previous works, and we identify the source of this error 
to be the presence of a finite divergence. 

Conclusions. We provide a method of quantifying the effect of a finite divergence in numeri- 
cal fields, together with detailed diagnostics of its sources. We also compare the efiiciency of 
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two divergence-cleaning techniques. Tliese results are applicable to a broad range of numerical 
realizations of magnetic fields. 

Key words. Magnetic fields, IVIethods: numerical, Sun: surface magnetism. Sun: corona 



1. Introduction 



Many astrophysical phenomen a, such as stellar and so l ar jets, flares, and corona l mass ejections 



are driven magnetically (e.g.. 



Tajima & Shibata 



2002 



Schriiver & Zwaan 



2008 



and references 



therein). Magnetically dominated plasmas are systems where the long-range, magnetic interaction 



dominates other f orces, e.g. 



low-corona {e.g.. 



Priest 



plasma pressure and gravi tational forces. A typical example is the 



2003 



Golub & Pasacho^koogh . There, the amount of energy associated 



with the magnetic field is much larger than other energy sources, and the dynamics of the coronal 



configuration is determined by the evolution of its magnetic field (e.g.. 



Forbes 



solar flares, where large c urrents develop in relatively small volumes (e.g., 



2000). This includes 



2011 



Shibata & Magara 



Aulanier et al.ll2012h . and coronal mass ejections (CMEs), which are powerful expul sions of 



coronal material that change the loc al configuration of the magnetic field drastically (e.g.. 



200C; 



Amari et al. 


2003 


Fan 


2010) 



Forbes 



20101) . In the coronal plasma, the magnetic energy is therefore the 



prime energy reservoir that fuels the dynamical evolution of these events. 

However, not all the magnetic energy is available for conversion into other forms of energy. 
Without changing the field significantly at the boundaries of the considered volume, the energy 
that can be converted into kinetic and thermal energies is given by the free energy, i.e., by the 
difference between the total magnetic energy and the energy of the corresponding current-free 
(potential) field. This very general result is known as Thomson's theorem, and it is based on the 
decomposition of the field into the sum of a current-carrying and a potential part. It does not depend 
on the presence of other forces, and is valid at any instant in time. 

The separation in the potential and free energies of Thomson's theorem is especially relevant 
for systems like the low-coronal field, that have different evolution time scales, as follows. The time 
scale of the coronal potential field is determined by the underlying photosphere, which is an inertia- 
dominated plasma, unlike the corona. This implies that the magnetic field at the photosphere has 
an evolution time scale that is much longer than the coronal one and that it is relatively insensitive 
to coronal changes. Since the magnetic field at the photosphere largely determines the coronal 
field's current-free component, the latter also evolves on the long photospheric time scale. As a 
consequence of Thomson's theorem, relatively fast events, suc h as flares and CME s , can only be 
powe red by converting part of the magnetic free energy (e.g. 



Aulanier et al. 



201C; 



Karpen et al 



20121) . 

In other words, the magnetic free energy is a sufficient condition for triggering active events, 
and it is considere d in the forecast of eruptions in the space weather context (see, e.g. 



Forbes et al 



2006 



Chen 



201 11) . Therefore, in this and similar applications, an accurate estimation of the free 
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energy is paramount for understanding the observe d magnetic field dynamics and the maxi mum 



energy that can be released in a flare or in a CME (lEmslie et al. 
references therein). 



2012 : 



Aulanier et al. 



2013 



and 



On the other hand, the free energy only provides an upper limit to the energy available for 



figuration does not need to be potential (see, e.g., 



Berger 


1985 


TavloJl986 


Low 


2001 



flare (reconnected) loops are frequently observed t o be sheared after a flare/eruption (see e.g., 



Asai et al. 



2003 



Lin et al. 



2010; 



Savage et al 



2012 



and re ferences therein), a feature that is also 



reproduced in numerical simulations (lAulanier et al. 



20121). This is an indication that post-event 



configurations have finite free energy, and the actual energy removed by the event is given by the 
diff^erence between the free energy of the pre- and post-event configurations. An assessment of the 
true energy budget related to a flare/eruptive event requires an accurate and reliable estimation of 
the magnetic energy. 

Another motivation of this study is to address the occurrence of unphysical magnetic config- 
urations. This is the case in some nonlinear force-free field (NLFFF) extrapolations when non- 
preprocessed, observed vector magnetograms are used as boundary conditions. The most obvious 
evidence of the nonphysical nature of some solutions is when the energy of the extrapolated field 
is lower than the potential field en ergy. This hap pens, for instance, in some of the solutions given 



in Table 3 of 



Metcalf et al 



( l2008h and Table 1 of 



Schriiver et al. 



( l2008h for three of the con sidered 



20101) . More 



extrapolation methods, including one used in the present manuscript (I Valori et al. 
generally, for all methods, the estimated coronal energy depends on the manipulations performed 
on t he observed data prior to their use in the actual e xtrapolation. (This step is called preprocess- 



ing. 



Wiegelmann et al 



2006 



Fuhrmann et al. 



2007 



) A significant part of the energy difference 



can eventually result from the details of the undergone preprocessing. 

As a result, the understanding of basic physical processes in the solar atmosphere requires an 
accurate estimations of the magnetic free energy. On the other hand, coronal models like NLFFF 
extrapolations, have shown that such accurate estimations are not easily obtained. In such cases, 
Thomson's theorem can be exploited to address the accuracy of (free) energy estimations. The 
fundamental assumption in Thomson's theorem is that the magnetic field is solenoidal. Such a 
property is only approximately fulfilled in numerical simulations and, more generally, in magnetic 
fields that are discretized on a mesh. A quantitative estimation of the effects caused by nonvanishing 
field divergence is complicated by its nonlocal nature. 

The main aim of this article is to quantify the effect of the presence of a nonsolenoidal compo- 
nent on the energy of a discretized magnetic field. This is studied using six different test magnetic 
fields that are a sample of the typical and characteristic examples used in the context of coronal 
solar physics. In the first part of the article, the energy of each test field is decomposed and in- 
terpreted using an extension of Thomson's theorem that can be applied to nonsolenoidal fields. In 
the second part we study how the energy changes, starting from a solenoidal version of each test 
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field and adding a parametric divergent component. The method and results of this study are of 
interest when working with any discretization of magnetic fields, e.g., for 3D coronal magnetic 
field extrapolations, as well as for magneto-hydrodynamic (MHD) simulations. 

In Sect.|2]the Thomson theorem for the energy of a magnetic field is summarized. The extension 
to nonsolenoidal, discretized fields is presented in Sect.|3] Section|4]introduces the six discretized 
fields together with their corresponding solenoidal versions that are used as test cases for applying 
our analysis, whose results are given in Sect. |5] Possible source of errors in our analysis are sort 
out in Sect. |6] Then, in Sect. [7] we present the parametric study of the energy dependence on the 
amount of divergence added to solenoidal magnetic fields. An analysis specific to numerical fields 
obtained by NLFFF extrapolations of observed vector magnetograms is presented in Sect. |8] and 
conclusions are finally given in Sect. |9] 



2. Magnetic energy of solenoidal fields 

We first consider the decomposition of the magnetic energy for perfectly solenoidal fields. By 
decomposing the field B as the sum of a potential. Bp = V0, and a current carrying contribution, 
Bj, 



B - Bp + Bj , 

the total magnetic energy E in CGS-Gaussian units in a volume is given by 

E s — f dV 
Stt 

= £p + £j + -L r (0Bj)-dS-^ r 0(V-Bj)dT, (1) 
where 

d'V represents the boundary of 'V, dS = n dS' , and n is the external normal to the bounding surface. 
Two conditions are classically considered: 



[a] n ■ (B - Bp)|fl^ - 0, i.e., the potential field Bp is computed from the same distribution of normal 
field of B on the boundary of 'V. This condition implies that n ■ Bj|a^ = and the surface 
integral vanishes in Eq. ([TJ; 

[b] V • Bj = 0, in which case also the rightmost volume integral in Eq. ([TJ vanishes. 
If these two conditions hold, then 



E^Ep+ E,, 



(2) 
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and the energy of a magnetic field is bounded from below by the energy of the corresponding 
potential field that has the same distribution of the normal component on the boundary of the 
considered volume. When applied to discretized fields, the above result holds under the implicit 
assumption that fields are numerically well resolved, yielding, in particular, continuous derivatives. 

The mathematica l equivalent of Eq. dJJl is known as Thomson's (or Dirichlet's) theorem, see 
e. e. . iLawrencd (1 1 998h . 



To satisfy the above requirement j^j, the scalar potential 4>{x,y,z) is computed as the solution 
of the Laplace equation 

f = 

■ (3) 

[ {dcj,ldh)\g^ = (n • Mirv 

In practical applications, Eq. Q can be solved numerically using standard methods. In the ap- 
plications presented in this paper, the Poisson solver included in the Intel® Mathematical Kernel 
Library was used. 



3. Magnetic energy of nonsolenoidal fields 

In this section we provide expressions for evaluating errors in the energy that stem from an imper- 
fect fulfillment of the solenoidal property, as is the case for discretized magnetic fields. In deriving 
Eq. ([U the divergence theorem, i.e., 

fv udO^^ r u dS, (4) 

is used, which may not be fulfilled by the techniques employed in constructing the numerical 
representations of magnetic fields or in their analysis. Moreover, if the numerically computed 
potential field Bp and current-carrying field Bj have a finite divergence, additional contributions 
can appear in the corresponding energy terms, Ep and Ej. 

We, therefore, seek a formulation of Eq. dl} for applications to numerical, nonsolenoidal fields 
that includes all possible sources of errors separately, that satisfies the requirement and that 
includes only volume integrals (thus avoiding using the divergence theorem). To obtain that, we 
first introduce the method of computing the potential and current-carrying parts. 

3.1. Helmholtz decomposition of tine potential part of tine fieid 

The accuracy in the numerical solution of Eq. ^ is limited, which may result in a finite divergence 
of the potential field. To quantify its effect, we can write 

[ AC = V ■ Bp 

Bp=Bp,, + Vf, where 1 ' (5) 

[ {dadh)\av = 
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which separates in Bp the solenoidal part, Bp s s Bp - Bp ns from the nonsolenoidal one, Bp ns = 
V^. This is equivalent to adopting the Hebnholtz decomposition for the vector Bp, together with 
the choice that all the nonsolenoidal component of Bp is contained in . Finally, the boundary 
condition for ({x, y, z) in Eq. (jSj is chosen such that Bp s satisfies the same boundary condition as 
Bp; i.e., they both fulfill requirement [a]. 

In practical applications, we first solve Eq. (|3]l numerically to determine <p, then we compute 
Bp = V0, and finally we overwrite the values of the normal components of Bp on each boundary 
according to Eq. (O. Since the latter operation enforces the requirement |ja], then any residual 
inaccuracy in the solution of Eq. Q, close to the boundary, implies a jump in the field, i.e., a finite 
divergence that adds to the divergence of the potential field discussed above. Second, we solve 
Eq. ^ to compute the residual nonsolenoidal component in Bp. 



3.2. Helmholtz decomposition of tine current-carrying part of tine fieid 

Using the Helmholtz decomposition on Bj we define a solenoidal component, Bj s, and a non- 
solenoidal one, Bj ns, such that 

i Ai/r = V ■ Bj 

B, =Bj,s + ViA, where { , (6) 

[ {di/f/dh)\s^ - 

the nonsolenoidal part of Bj being: Bj ns = Vi/'. The boundary condition for ifr in Eq. ^ is chosen 
to have the same boundary condition for Bj s and Bj, i.e., to fulfill the requirement |jaj. As for the 
potential field, the required values of Bj ^ at the boundaries (i.e., zero in this case) are overwritten 
onto the solution of Eq. (|6]l which is obtained numerically, so that any error in matching these 
values by t//{x, y, z) reduces to a finite jump close to the boundaries. 

Finally, we notice tha t this method is often used to remove the divergence of a vector field 



jBrackbill & Barnesll 19801 sometimes referred to as "projection method"), and it has the property 



of conserving the current, i.e., V x Bj = V x Bj ^. 



3.3. Gauge-invariant decomposition of tine magnetic energy 

We now summarize the procedure for the decomposition of the magnetic field. For a given nu- 
merical magnetic field B, we solve Eq. (O numerically and compute the corresponding potential 
component Bp and current-carrying component Bj - B - Bp. Next, we compute the solenoidal 
component Bp s = Bp - and the nonsolenoidal component Bp ns - of the potential field by 
solving Eq. (|5]l numerically. Similarly, the numerical solution of Eq. (|6]l provides the solenoidal 
component Bj.s = Bj - Vi/* and the nonsolenoidal component Bj ^s - Vf/r of the current-carrying 
part of B. The values of the different components at the boundary are such that the condition [a] 
is satisfied (Sect. |2]i. Finally, by substituting the field decomposition in E = ]^ B^jJ^/Stt and 
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grouping it again as in Eq. ([T]i, we obtain 



with 



j^Vi-Wifrd^+ Bp,, ■ Bj,s d^j . 



V J'V 

^ (8) 



All terms in Eq. O are positively defined, except for Zsmix- For a perfectly solenoidal field, it is 
Ep^s = Ep, Ej.s = Ej, Ep^ns = i'j.ns = ^mix = 0, and Eq. Q reduces to Eq. ©. 

Finally, Eq. is normalized such that 

1 ~ ^p,s "t" ^J,s "t" ^p,ns "t" ^J,ns "t" ^mix 5 (9) 

where the tilde indicates that the corresponding definition in Eq. ((SJ is divided by E. 

Using the divergence theorem, Eq. (|4]i, and the condition several terms in the above expres- 
sions could be simplified. However, since practical test fields may be obtained with methods that 
do not insure that the divergence theorem holds numerically, we have kept all the terms in Eq. 
Indeed, the simplification obtained by using the divergence theorem results in mixing other numer- 
ical issues with the issue of the finite divergence, producing cumbersome results, up to the point 
where Eq. (|7]i is not satisfied numerically. Moreover, the direct appearance in the integrals of the 
scalar potentials, rather then their gradients, introduces an undesired gauge-dependence. 



3.4. Sources of the violation of tiie Thomson theorem 

We summarize which are the source of errors that we consider in Eq. Q. First, the energy is 
aff'ected by the finite divergence of the current-carrying part of the magnetic field, which enters the 
Zij ns and ^niix terms. Additionally, the potential field may have a finite divergence, owing to the 
Hmited numerical accuracy of the solution of Eq. (|3]l, both in the volume and close to its boundary. 
These effects are contained in the Ep^^s and E^ix terms. 

As long as these are the only source of errors, then the sum of the terms on the righthand side 
of Eq. (I2I) must be equal to the total energy E computed using B directly, and Eq. (|7]l must hold 
numerically even for nondivergence-free fields. Equivalently but using normalized quantities, the 
sum on the righthand side of Eq. (|9j must be equal to one. We show in Sect.|5]that the total energy 
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is indeed retrieved by the decomposition we adopted, allowing us to identify the source and extent 
of the eventual violation of Thomson's theorem, Eq. (|2]i. 

3.5. Accuracy of the decomposition of tine energy equation 

A further step is the assessment of the accuracy of the decomposition, Eq. (|7]). First, we address 
how effective the decomposition in the solenoidal and nonsolenoidal parts is in concrete numerical 
applications in Sect. 16.11 

Second, the continuity condition, implicit in the derivation of Eq. Q, implies that numerical 
derivatives can be computed precisely enough in the employed discretization. This may not be the 
case in some numerical applications, e.g., when observed values are used as boundary conditions 
for computing magnetic fields. The continuity of the fields in relation to small scales is discussed 
in Sect.[8l 

Finally, our decomposition employs the numerical solution of Laplace and Poisson equations. 
We briefly recall the conditions for uniqueness of the general Poisson equation 



[ idu/dh)\g^ = g 

where f{x, y, z) is a source term in "V, and g is the boundary value on d'V. The use of Neumann 
boundary conditions imphes that the solution u{x, y, z) is only unique up to an additive constant. 
For Eqs. (O, (|5]l and the freedom in the additive constant is equivalent to a gauge freedom 
for the scalar potentials (p, and i/r, respectively. This gauge dependence is, however, irrelevant 
for Eq. O, since the energy decomposition is intentionally derived in a way such that the scalar 
potentials only appear in conjunction with the gradient operator. 

Integrating Eq. (fTOl i in 'V and using the divergence theorem, Eq. (|4]i, we find that source and 
boundary values must satisfy 



which is a necessary condition for the uniqueness of the solution u. This implies that, for Eq. (O 
where f - and g - n - Rlyy, the flux of B through d'V must vanish. For Eq. (|5]l, where / = V ■ Bp 
and g - 0, it implies that the volume integral of V ■ Bp must vanish. Similarly, for Eq. (|6}, where 
/ = V • Bj and g - 0, the volume integral of V ■ Bj must vanish. When such conditions cannot be 
insured, the uniqueness of the solution is not guaranteed. The effect of the violation of Eq. (fTTT i is 
studied in Sect. 16.21 

4. Test fields 

To explore the effects of a finite divergence in a representative sample of practical situations, we 
consider six test fields Btest obtained from analytical models, numerical simulations, and NLFFF 




(10) 




(11) 
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Fig. 1. Selected field lines of the six test cases: (a) the potential field of a double dipole, Bdd; (b) the 
TD model, Btd; (c) the MHD model, Bmhd; (d) the NLFFF model of the nonpreprocessed magnetogram of 
AR 1 1 158, Bexi ; (e) the NLFFF model of the preprocessed magnetogram of AR 1 1024, Bex2pp; (f) the NLFFF 
model of the nonpreprocessed magnetogram of AR 1 1024, B^xi- The vertical component of the magnetic field 
at the bottom boundaiy is shown on a gray scale, with the positive (respectively, negative) polarity in white 
(respectively, black). The different line colors outline different types of connectivities. 



extrapolations. Their magnetic configuration is outlined in the field-line plots in Fig. [T| Further- 
more, we consider six additional test cases Btest,s, which are obtained from each of the Btest by 
removing the nonsolenoidal part of the field. 



4. 1. Discretized analytical test fields 



The first test field that we consider is the potential field Btest = Bp p generated by a pair of vertical 



Torok & KliemI (l2003h for the 



magnetic dipoles, located at iO,+yDD,ZDD), see, e.g., Eq. (7) in 
analytical expression of the field. We set yoo - 2 and zdd = -1.5, and the field is normalized such 
that the z-component has a maximum value equal to unity at the bottom boundary (z = 0). The 
only currents and finite divergence errors present in Bdd are generated by truncation errors in its 
discretization. 

The second e mployed test field, Btest = Btd, is the model of the magnetic field of an active 



region derived in 



Titov & DemoulinI (Il999h . given by a section of a curre nt ring surrounded b y a 



Valori et al 



(2012), to 



stabilizing potential field. The employed configuration is the same as in 
which we refer the reader for further details. In this case, the test field has an explicit current- 
carrying component sustained by a flux rope. The analytical formulae defining the test field are 
approximate, which together with the rather coarse resolution employed here, yield relatively large 
finite-divergence errors. 

For bofli test fields Bdd and Btd the discretized volume is = [-12, 12] x [-19, 19] x [0, 16], 
with uniform resolution A = 0.12 in all directions. 



4.2. Numerical tests fields 



The next test field that we consider, Btest = Bmhd, is a snapshot of a m agneto-hydrodynami c nu- 
merical simulation of magnetic reconnection in a null-point topology (iMasson et alj|2012l) . To 
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use our present-stage diagnostic, we interpolated the original snapshot onto a uniform and homo- 
geneous grid, whereas the original simulation was performed using a nonuniform one. Because 
the divergence values are slightly incr eased by the interpol ation, they are not representative of the 



quaUty of the simulations presented in 



Masson et al. 



( l2012h . However, they still serve our purpose 



of providing a typical situation arising from the numerical evolution of magneto-hydrodynamic 
equations. The considered volume is = [-20,0] x [-20, 10] x [0, 12] with uniform resolution 
A = 0.05 in all directions, and the field is normalized such that the vertical component is unity at 
its maximum. 

Next, we consider three NLFFF ext rapolations o f Hino de/SOT vector magnetograms, obtained 



with the magneto-frictional method in 



Valori et al 



(12010 ). The original resolution of the vector 



magnetograms is 0.3", and they can be preprocessed (iFuhrmann et al.ll201 11) to improve their com- 
patibility with the force-free assumption on which the extrapolation code is based. 

Our fourth test field, Btest = BexI, is the nonlinear extrapolation of a vector magnetogram of 
AR 1 1 158, measured on 14 February 201 1. The vector magnetogram was binned to the resolution 
A = 1.1" prior to extrapolation, and no preprocessing was applied in this case. The analyzed 
coronal model volume in arcsec is ^ = [-21, 68] x [-273, -171] x [0, 123]. The Hinode/SOT field 
of view of the measurements employed for this extrapolation cuts through the external sunspots of a 
quadrupolar field distribution, resulting in high field values at the lateral edges of the magnetogram. 
Even computing the potential field is problematic in this case, therefore we limited the considered 
volume to the bipolar core of the extrapolated field. 

The fifth test field, Btest = Bex2pp, is the extrapolated coronal field model above AR 1 1024 on 
4 July 2009. In this case, the full resolution of Hinode/SOT is used, and the vector magnetogram 
is preprocessed before extrapolation. The extrapolation covers a volume of ^ = [-41,42] x 
[-141, -16] X [0, 98] arcsec, with uniform reso l ution A - 0.3". This model of the coronal field of 
AR 11024 is discussed in detail in 



Valori et al 



(1201 lb . where more details about extrapolation of 



vector magnetograms can be found. 

Finally, the sixth test field, Btest = Bex2, is the same case as Bex2pp except that the vector mag- 
netogram is not preprocessed prior to extrapolation. More details on the numerical implementation 
are given in Appendix lAl 



4.3. Cleaned test fields 

Since a small divergence of B is one major condition for the Thomson theorem, Eq. ([2}, for each 
test field Btest we consider a corresponding solenoidal version of it, Btest.s, which is computed from 
Btest employing the divergence cleaner described in Appendix iBl In Cartesian coordinates, such a 
solenoidal field has the same x- and y-components as Btest, whereas the z-component is changed 
everywhere in the volume, except for the top boundary. Therefore, Btest and Btest.s have the same 
distribution of normal field on all boundaries except for the bottom one, where Btest.s differs from 
Btest by an amount that is related to the combined effect of V • Btest in the whole volume. Since the 
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Table 1. Numerical tests of Tliomson's theorem. 



■ 1 

i»test 


mm} 


E 




£'j,s 


^p,ns 




£'j,ns 








Sum 


■1 

»DD 


2x 10 


-6 


1 /I ^ 


i.UU 


U.UU 


4 X iU 


-5 


/ X iU 


-6 


—4 X iU 


-3 


1 r\r\ 
i.UU 


Btd 


3 X 10 


-6 


3.90 


0.81 


0.16 


3 X 10 


-5 


6x 10- 


-4 


0.02 




i .UU 


Bmhd 


2x 10 


-5 


1.94 


0.94 


0.06 


1 X 10 


-6 


1 X 10" 


-4 


-1 X 10 


-3 


1 MM 

i.UU 


BexI 


4x 10 


-3 


4.21 


0.79 


0.38 


5 X 10 


-4 


0.29 




-0.46 




1 MM 
i.UU 


Bex2PP 


9x 10 


-4 


1.51 


0.88 


0.11 


2x 10 


-4 


0.14 




-0.12 




1 MM 
i.UU 


Bex2 


2x 10 


-3 


0.72 


2.29 


0.14 


3 X 10 


-4 


0.94 




-2.38 




1 MM 
i.UU 


Btest.s 






Bdd.s 


1 X 10" 


IS 


1.44 


1.00 


0.00 


4x 10 


-5 


5 X 10- 


-5 


7x 10- 


4 


1 MM 
i .UU 


Btd.s 


4x 10- 


21 


3.95 


0.84 


0.16 


3 X 10 


-5 


4x 10- 


-5 


-3 X 10 


-4 


1.00 


Bmhd.s 


3 X 10- 


21 


1.94 


0.94 


0.06 


1 X 10 


-6 


1 X 10- 


-6 


-1 X 10 


-5 


1.00 


BexI.s 


6x 10- 


18 


5.98 


0.43 


0.57 


2x 10 


-4 


5 X 10- 


-3 


-6 X 10 


-3 


1.00 


Bex2PP,s 


2x 10- 


17 


3.15 


0.42 


0.58 


1 X 10 


-4 


2x 10- 


-3 


-3 X 10 


-3 


1.00 


Bex2,s 


8 X 10- 


18 


0.99 


0.61 


0.39 


2x 10" 


-4 


1 X 10- 


-3 


-2 X 10- 


-3 


1.00 



Notes. The employed test fields, defined in Sect.|4l are named in the leftmost column. Second column, ( |/| ): 
the divergence metric of the fields (see Eq. iC.H ). Third column, E: energy of the test fields in units of 
10^^ erg. The Bdd> Btd, Bmhd fields (and their coiTesponding solenoidal fields Bdd.s. Btd.s, Bmhd.s) were 
rescaled assuming a maximum value of the photospheric vertical field equal to 300 G and a typical distance 
between the sunspot's centers of (50, 50, 120) Mm, respectively. The successive five columns are the difi"erent 
contributions to Eq. (|9}, and "Sum" corresponds to their sum. All terms from "£p s " to "Sum" are normalized 
by E. £p_s is the magnetic energy of the potential field Bp_s, £j.s that of the solenoidal component of the current- 
carrying one Bj_s, £p „s and £j „s are the contributions associated to the divergence of Bp and Bj, respectively, 
£,„ix is a mixed potential-current carrying term (see Eq. (|8) for their definitions). 



divergence cleaner changes the value of the normal field component on one boundary, the potential 
fields computed from the boundary values of Btest and of the corresponding solenoidal Btest.s are 
not the same. Additionally, the divergence cleaner alters the current of the field, as prescribed 
by Eq. (IB.5I ), of an amount that is proportional to the divergence of Btest- Therefore, the field 
that is obtained by applying the cleaner may have drastically different properties than the original 
field. Finally, let us notice that different solenoidal fields can be derived from Btest using different 
methods. The divergence-cleaned versions of the test fields Btest.s are used here as illustrative 
examples. 



5. Numerical tests of Thomson's theorem 

In this section we apply Eq. (|9| to the test cases described in Sect.|4] Table[T]summarizes the values 
of the divergence metric ( | ) defined in AppendixICland the contribution of each term to Eq. (|9]l, 
for all test fields. The divergence metric spans values from 10-^' to lO-"'. In all cases, the rightmost 
column, corresponding to the sum of the righthand side of Eq. (|9]l, is equal to unity, despite the large 
difference in the divergence values. Therefore, we conclude that Eq. (|9]i completely accounts for all 
relevant contributions to the energy, in all test cases. We then consider the different contributions 
to Eq. (|9]i case by case. 
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5.1. Results with the test fields 

The top part of Table [T]refers to the test fields Btest- In general, the energies E of the different test 
fields go from the purely potential case of Bdd, where s = 1, to high-free-energy cases (Btd 
and Bexi, with Ep^^ - 0.8), where the field is strongly nonpotential. The main source of violation 
of Thomson's theorem, Eq. (|2|i, in all cases is the mixed current-potential term E„iix, except for the 
Bex2pp case where £j ns is slightly higher in absolute value than £niix- 

More precisely, Bdd is nearly perfectly potential, with nonsolenoidal spurious fluctuations con- 
tributing to the total energy for few parts per thousand at most (in £mix)- Bjd has a 16%-energy 
contribution from the current-carrying part of the field Ej^^, with a 2% contribution from the non- 
solenoidal field related to the current-carrying structure (in E^i^ but not in £j,ns)- This is the effect 
of the approximate nature in the matching between current-carrying and external potential fields in 
the equilibrium that defines the Btd field. Bmhd, which has 6% free energy £j s, has an even lower 
nonsolenoidal contribution (-0.1%). In all three cases, there is very small (Btd) or no significant 
(Bdd, Bmhd) violation of Thomson's theorem. 

We now move to the NLFFF extrapolations. These show values of ( |//| ), which are two-to-three 
orders of magnitude greater than in the first three cases. The contribution of the nonsolenoidal part 
of the potential field to the total energy, Ep^„^, is always negligible with respect to the other terms. In 
the Bex2pp case, the free energy associated with the solenoidal part of the current-carrying field ^ 
is about 11%, and the potential field energy is 88% of the total energy. The sum of the potential 
and current-carrying solenoidal parts accounts for 99% of the total energy, apparently verifying 
Thomson's theorem accurately. However, £j ns is 14% and fimix is -12%; i.e., the errors related to 
the divergence of the current-carrying part of the field have comparable magnitudes and compensate 
for each other. These are the dominant sources of error, almost three orders of magnitude more than 

The test case with the highest value of ( \ fi\ ) is Bexi. With respect to the Bex2pp case, BexI is 
characterized by three times higher free energy £j_s, twice the error on current Es „s, and almost 
a four times larger error on ^mix. Again, the last two are largely compensating each other. We 
conclude that the interpolation to one third of the resolution used for BexI is less efficient than 
preprocessing (used for BexIPp) in eliminating the source of violation of Thomson's theorem. 

This situation is even more extreme in the case of the extrapolation of the nonpreprocessed, 
noninterpolated magnetogram Bex2- Although this case has a value of the mean divergence {\fi\} 
that is only a factor two higher than for Bex2pp, and not even the highest one, it shows the most 
pathological behavior: The potential field has an energy 2.29 times the energy of the test field, 
which is downright unphysical according to Eq. ([2}. Such a high value is compensated for by an 
equally high value of E,^a (-2.38). On the other hand, the current-carrying part of the field Ej^^ 
accounts for 14% of the energy, but the associated error Ej.ns is more than six times larger. Such 
large errors are related to the high values of the divergence — in particular at the bottom boundary — 
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and their actual values are very sensitive to the numerical details of the computation. Additional 
analysis of Bex2 and Bex2pp is discussed in Sect. 18.11 

We finally notice that in the preprocessed case Bex2pp, the error from £j_ns or might be 
considered as still tolerable if compared with the total energy (errors on vector magnetograms are 
similar, after all), but it seriously compromises the reliabihty of the free energy estimation, each 
one being as high as fij s itself. 

5.2. Results with the cleaned test fields 

We now consider the bottom part of Table [T]for the solenoidal fields. The values of the divergence 
are drastically reduced in all cases to 10 or less, which shows that the cleaner in Appendix iBlis 
an eff'ective — and fast — way of removing the nonsolenoidal component of a discretized magnetic 
field. For the purpose of this article, we can then consider all Btest.s to numerically be perfectly 
solenoidal. All error terms, i.e., Ep ns, £j,ns, and E,^a, are smaller than 1%, and we recover Eq. ([2) 
in a numerical sense. 

More precisely, the Bqd.s and Bmhd.s cases are practically identical to their corresponding test 
fields, as far as the energy metrics E, Ep.s, and £j s are concerned. On the other hand, Btd.s shows 
an increase of about 1.3% of the total energy, E, as a result of the removal of the error in E^a 
of Btd- The error removal affects the potential field energy Ep more, which raises about 5% with 
respect to the energy of the potential field in Btd (in nonnormalized values), as a consequence 
of the cleaner's modification of the bottom boundary. In contrast, the relative contribution of the 
current-carrying part £] s is unaffected by the cleaner. It is true that ( \fi\ ) differs by 15 orders 
of magnitude between Bjd and Btd.s, but it is significant anyway that the removal of a 2%-error 
in fimix changes the nonnormalized values of the total energy E and potential field energy /Sp ^ 
of 1% and 5%, respectively. We conclude that, even in relatively divergence-free fields, residual 
nonsolenoidal efifects can be energetically significant. 

In the extrapolated cases, the removal of the larger divergence has far stronger consequences. 
In the first place, the nonnormalized field energy E of the cleaned fields Bexi.s, Bex2pp,s, and Bex2,s 
is increased of 42%, 109%, and 38%, respectively, with respect to those of the corresponding test 
fields. As a consequence of the higher values of E, the importance of potential fields relative to 
the total energy ^ is decreased (to 78%, 95%, and 40% of their test-field values, respectively). 
In contrast, the energy contribution related to the current-carrying part of the field ^ is strongly 
increased, as expected, since the cleaner introduces currents that are related to the cumulated di- 
vergence that is removed (see Appendix IbTi. 

We conclude that the cleaned fields that are obtained from the test ones using the method in 
Appendix |B] all comply with Thomson's theorem accurately. However, three of them, namely 
BexI.s, Bex2pp,s, and Bex2,s, are energetically very different from the original fields Bexi, Bex2pp, 
and Bex2, respectively. Incidentally, we notice that the removal of the finite divergence does not 
conserve the approximate force-freeness of the extrapolated fields. 
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Table 2. Contributions to f^ix in Eqs. I |8l9t . 



Btest 


-^p,s/p,ns 


E],s/J,ns 


■^p,s/J,ns 


■^J,s/p,ns 


■^p,ns/J,ns 


^P,s/J,s 


F ■ 


Bdd 


-0.01 


0.00 


0.00 


-0.00 


-0.00 


0.00 


-0.00 


Btd 


-0.01 


0.00 


0.03 


0.00 


-0.00 


-0.00 


0.02 


Bmhd 


-0.00 


0.00 


-0.00 


-0.00 


-0.00 


-0.00 


-0.00 


BexI 


-0.01 


-0.08 


-0.43 


0.00 


0.00 


0.05 


-0.46 


Bex2PP 


0.00 


-0.03 


-0.10 


0.00 


0.00 


0.01 


-0.12 


Bex2 


0.00 


-0.28 


-2.46 


0.00 


-0.00 


0.36 


-2.38 



Notes. £p.,/p,„, = 1 Bp,, • d'V, E,,p,,, = ^ Bj,, • d'V, 
£p,s/j,„s = ^j^ Bp, ■ d'V, £j,/p.„, = i Bj,, • d^, 

4.ns/J,n,s = I ■ d^, £p,,/j,, = i Bp, • Bj., d^, 

5.3. Contributions to E^ix for tiie test fields 

In many of the test fields in Table [1] E^i^ is the largest source of error. Table shows the six 
contributions to Emi^ in the order in which they appear in Eq. ([8} and their sum Emi^ for the six 
test cases Btest. We do not consider the solenoidal fields Btest.s since all terms are mostly zero and 
never bigger than 0.7% The following conclusions can be drawn. First, the main contribution to 
£mix is £p,s/j,ns = ■g iy ^p,s ' ^lA d'V in all cases. The main source of violation of Thomson's 
theorem is then the divergence of the current-carrying part of the field. More often than not, this 
term has a similar magnitude and opposite sign of Ej ^^, which is positive-definite. However, there 
is no obvious reason for £p,s/j,ns to be always — or predominantly-negative, and we regard this as a 
coincidence. 

Second, the terms with residual divergence of the potential field (i.e., any term containing 
in Eq. ([8]l) are always negligible. Therefore, also in view of the always low £p ns values in Table [1] 
we can conclude that the divergence of the potential field always gives a negligible contribution to 
the energy. 

Third, the integral Ej.s/j.ns = jy Bj„s ■ Vi/r d'V, and the integral -Bp.s/j.s = jy Bp_s ■ Bj,s d'V 
have finite values for the extrapolations in Table |2] Analytically, they should be vanishing. Using 
the divergence theorem, Eq. (|4), the surface integral vanishes because Bj sl^)^ = 0, and the volume 
integral vanishes because V-Bj,s - 0. The first condition is enforced at the boundary, but the second 
is only approximately true numerically (see also Sect. l6.lt . This is not enough to insure that fij.s/j.ns 
and fip.s/j.s vanish numerically. This is why we adopted the decomposition of the energy of Sect. 13. 31 
that only contains volume integrals. 

6. Source of errors in the decomposition 

6.1. Values of { \fi\ ) for the field decomposition in Eqs. (3ltBi 

In this section we quantify how accurate the decomposition in solenoidal and nonsolenoidal con- 
tributions is. Table [3] reports the values of the logarithm of ( |/| ), defined by Eq. ( IC.2b . for the 
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Table 3. Values of logjQ(( |/| », for the fields decomposition in Eqs. ([3}, (|5}. and ((SJl (see Eq. iC.2i for ttie 
definition of ( |/|| )). 



R. , 




R 


R 

"p.s 


R 




Rt 


Ri 


Bdd 


-5.61 


-4.98 


-5.60 


-2.43 


-2.35 


-2.93 


-2.46 


Btd 


-5.54 


-4.84 


-5.40 


-2.29 


-4.22 


-4.71 


-3.20 


Bmhd 


-4.78 


-5.76 


-6.23 


-2.55 


-3.96 


-4.54 


-2.45 


BexI 


-2.40 


-2.60 


-2.66 


-1.42 


-2.08 


-2.36 


-2.04 


Bex2PP 


-3.05 


-4.02 


-4.09 


-2.28 


-2.62 


-2.90 


-2.62 


Bex2 


-2.66 


-3.87 


-3.96 


-2.12 


-2.69 


-2.86 


-2.80 


Btest.s 


Bdd.s 


-18.0 


-4.98 


-5.60 


-2.44 


-0.83 


-2.62 


-1.91 


Btd.s 


-20.4 


-4.84 


-5.40 


-2.31 


-3.41 


-4.09 


-2.15 


Bmhd.s 


-20.5 


-5.42 


-5.79 


-2.39 


-4.39 


-4.70 


-2.31 


BexI.s 


-17.2 


-2.65 


-2.71 


-1.47 


-1.72 


-2.11 


-1.83 


Bex2PP,s 


-16.8 


-3.78 


-3.93 


-2.26 


-0.44 


-2.76 


-2.25 


Bex2,s 


-17.1 


-3.63 


-3.79 


-1.99 


-1.51 


-2.71 


-2.19 



Notes. Column B here is the logarithm of the column { |/(B)| > in Table [T] More negative values correspond 
to more solenoidal fields. 

field decomposition used in Eq. O. Since ( |//| ) is not additive in the field, its value for, say, B is 
different from the sum of its values for the potential Bp and current-carrying Bj components. 

We next consider the decomposition of the potential field given by Eq. ^ for the test fields 
Btest (upper half of Table |3]l. Values of ( | ) for the solenoidal part of the potential field Bp 5 are 
better {i.e., more negative) than those for Bp ns, so that the Bp s is indeed more solenoidal than Bp ns. 
However, it is only in the first three cases, Bdd, Btd, and Bmhd, that logio(( \fi\ )) has a noticeably 
more negative value for Bp ^ than for Bp. In the other cases, the values are relatively close to each 
other, and Bp s is only marginally more solenoidal than Bp. On the other hand, Bp ns is always much 
less solenoidal than Bp. This is partly the effect of the nonadditivity of the metric { \fi\ ), and partly 
because Bp ns is, on average, much smaller than Bp s, as the corresponding energy metrics in Table [T| 
show. (In particular, fip.ns, i-e., the energy associated with Bp ^s, is always extremely small.) 

Similar conclusions can be drawn looking at the decomposition of the current-carrying part, 
Bj, where this time the energy associated with the solenoidal error (see Ej ns in Table [T]) is more 
significant. In this case, values of ( ) for all three contributions Bj, Bj s, and Bj „s are of similar 
magnitude. Again, the nonsolenoidal part, Bj ^s, has a higher divergence value than the solenoidal 
one, Bj s, but only marginally so for Bdd and extrapolated fields. 

We consider the solenoidal test fields Btest.s (bottom half of TableO. Values of ( | ) for a given 
field component belonging to Btest and to the corresponding Btest.s are very similar. For instance, 
the value of log[Q(( \ fi\ )) for Bj s in, say, the test field Bexi is -2.36, whereas for the corresponding 
contribution for Bexi,s it is -2.11. Therefore, the above discussion of the contributions to Btest 
holds for those of Btest.s as well. In contrast, the total divergence of the field is very different in 
the two cases, i.e., -2.4 and -17.2, respectively. This is a clear indication that the accuracy of the 
field decomposition is determined by the accuracy in the solution of Eqs. (l3]|6l) rather than by the 
divergence of the total field. 
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Table 4. Removal of flux imbalance. 



B, 



'test 



<|/;(B[g,)|) OaT.(Btest) 1>a^(Blg.) £q> 



Bdd 

Btd 
Bmhd 

BexI 



-5.61 -6.83 -26.4 -17.3 -11.2 

-5.54 -7.79 -26.0 -19.1 -11.8 

-4.78 -4.36 -22.8 -12.2 -6.29 

-2.40 -1.33 -18.7 -5.18 -3.30 

-2.60 -2.00 -20.9 -7.49 -4.52 

-2.63 -1.86 -20.0 -6.77 -4.35 




Notes. For all quantities, the logjg of the absolute value is shown. ( |/;| > is the divergence metric defined in 
Eq. l lC.2t : <i>g^> is the normalized magnetic flux through all boundaries, Eq. il2i : Ef, and £o,mix are normalized 
to the energy E of Bl'^Jp which is the flux balanced field associated to Biest, see Eq. l ll3t . 

In conclusion, the Poisson solver provides a decomposition of the magnetic field where the 
solenoidal parts have a smaller divergence than the original field, as required. The limit in the 
accuracy of the decomposition comes from the accuracy of the solver, and not from the level of 
solenoidality of the initial field. One possible source of inaccuracy for the solver is the incompati- 
bility of the boundary conditions used in Eqs. (|3]|6]l, which is discussed in the next section. 

6.2. Compatibility of boundary conditions in Eqs. (3l6i . 

We here consider the normalized flux of the field, 0,9^, computed as the surface flux through all 
six boundaries, normalized to the mean flux entering and exiting from the lower boundary: 



The values of logjQ |4>(5'v(Btest)l in Table|4]show that the test fields of the extrapolation cases Bexi, 
Bex2pp, and Bex2 are not flux-balanced. Therefore, the decomposition of Eq. ^ based on the 
solutions of Eqs. (|3]|6]l may be inconsistent (see Eq. (fTTT i and related text). The purpose here is to 
determine whether the unbalanced flux affects the accuracy of any of the terms in Eq. (|7]i. 

A flux-balanced field, R^^^^, can be computed from a flux-unbalanced one, Btest, by splitting the 
original field as 



and assuming Bo = V0 to be generated by an uniformly distributed, constant divergence; i.e., 
A0 - constant. We choose the simple solution oc r^, and fix the constant such that the flux of 
Bo equals the flux of Btest, yielding 



Table |4] shows that the above modifications to Btest is effective, drastically reducing the net 
flux of the original field to very low values (compare logjg |'tc)'v(Btest)l with logn, |<l>f}^(B]i|J!,\)|). On 
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(12) 



Bte^st - Btest + Bo , 



(13) 
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the other hand, the effect on the field of Bo is very small. Both energy terms related to that (i.e., 
Eq) = 4e ^o'^'^ and ^o.mix = ^ JyBtesfBod'V) are negligible (with a contribution below 0.01% 
of the total energy E computed for Bj^.^^). 

Repeating the same analysis of Sects. |5] and W\\ for the flux-balanced part of the field only, 
J^test' yields no significant change: all values in Tables[T]|2] and[3]are identical. Inaccuracies of the 
Poisson solver in solving Eq. (|3]l are therefore related to the solver itself, not to the incompatibility 
of the boundary conditions. 

In a similar way, the test field can be modified to have vanishing volume divergence, which is 
the requirement for consistency in solving Eqs. (1516) . using 



The result is likewise clear: no significant change is found in the values of Tables [T] |2] and|3] 

Therefore, an imperfect consistency of source and boundary conditions play no role in the 
accuracy of the solution of the Laplace and Poisson equations employed in the decomposition, 
Eq. (|7]i, for any of the test cases. Recalling the results of Sect. 16.11 we conclude that the accuracy 
limitation of our analysis comes from the solver itself. In this respect, we note that, when the 
method used in Eqs. (|5]l and ^ is viewed as an algorithm for removing the divergence (Projection 
method), it is far less efficient than our divergence cleaner described in Appendix iBl On the other 
hand, the projection method has other advantages; for instance, it change neither the current in the 
volume nor the normal component of the original field at the boundaries. 



7. Parametric study 

In this section we study how the relative energy of the field depends on its divergence in pro- 
gressively going from a solenoidal to nonsolenoidal realizations. The purpose is to offer a practical 
method of fixing the level of solenoidal eiTors that can be tolerated in a given numerical realization, 
based on their consequences on the energy of the field. 

7. 1 . Parametric models of finite-divergence fields. 

For a given test magnetic field Btest, the corresponding solenoidal field Btest.s is considered. A 
parametric, nonsolenoidal field B^ is obtained by adding a nonsolenoidal component Bjiv to Btest.s, 
using a control parameter 6, as 




B;; - Btest.s + ^ ^div ■ 



(14) 
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We consider here two models of Bdiv, namely 



BhIv — 



-z J~\W ■ BtestW Modell, 
4(xf^dy + yf d/+ (15) 
+z J^- dz') V ■ Btest Model 2. 



Adding the first divergence model for 5 = 1 is the inverse operation of the cleaner in Section |B] 
since B^(5 = 1) = Btest- For other 6 values, the resulting field B^ only differs from the solenoidal 
field Btest.s in the z-component. The second model for the divergence is more general, because it 
changes all three components of Btest.s in the volume, although not on all boundaries. 

Both divergence models in Eq. (fTST i are based on the computed V ■ Btest- In this way, we 
relate the divergence models to the source of error that is specific to the considered test field. For 
instance, we expect that errors in the test case Bdd, which are only generated by truncation errors, 
have a different distribution in space than those coming from the approximate nature of the Btd 
equilibrium, or from a numerically constructed field like Bex2pp- 

The influence of a finite divergence of B^ on the energy value can be written as 

E = ^^^div + 2(5£'s,div + -EtesLs , (16) 

where E, E^n, and /itest.s are defined as usual as proportional to the volume integrals of B^, B^^^, 
and B~^^^ ,., respectively, and Es.div = jy(Btest,s ■ Bdiv) dT/8;r. 

Below, the energy dependence on 6 is studied for the test fields described in Sect. |4] Since 
the separation in solenoidal and nonsolenoidal components is known by construction, we simplify 
the presentation by analyzing the energies of the total fields according to Eq. (fT6l l. and we do not 
separate the error sources as in Sect. |5] For each value of 6, we consider B^ as the test field to 
analyze and compute the corresponding potential field according to Eq. (O. 



7.2. Parametric dependence of the energy 

Figure |2] shows the energy for the two divergence models in Eq. ( fTsT i, as a function of the control 
parameter 5 in a wide range of values. Due to the large difference in V-Btest between the six models, 
the top- and bottom-rows have different scales. The orange lines show the energy normalized 
with the energy of Btest, which is not dependent on 6: They follow the expected parabolic profile 
of Eq. (fT6b . only scaled by the normalization factor. Model 1 (continuous orange lines) yields 
a smaller variation of the energy with 6 (corresponding to lower values of £div) with respect to 
Model 2 (dashed orange lines), and is centered farther away from 6-0 [i.e., Model 1 has higher 
values of E's.div/^'div)- 

The orange curves in the top row of Figure |2]show that it takes very high values of 5 in order to 
have a variation of order one of the energy in the Bdd, Btd, and Bmhd cases [e.g., for the Model 2 
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c) Bmhd 
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Fig. 2. Magnetic energy, Eq. l ll6t , normalized to tlie energy of the potential field having the same distribution 
of normal field on the boundaries (black lines) or to the energy of the reference field (orange lines), as a 
function of the amplitude S of the nonsolenoidal term, Eq. ( I14t . Two models for the divergence are shown 
according to Eq. dlSt . Each panel shows the results of a test field: (a) the potential field of a double dipole, 
Bdd; (b) the TD model, Btd; (c) the MHD model, Bmhd; (d) the NLFFF model of the nonpreprocessed 
magnetogram of AR 11158, Bexi; (e) the NLFFF model of the preprocessed magnetogram of AR 11024, 
Bex2pp; (f) the NLFFF model of the nonpreprocessed magnetogram of AR 11024, B^xi- An important change 
of scale of both axes is present between the top and the bottom rows. 



applied to Btd S-t S - 15). On the other hand, the energy of the extrapolated fields shows a much 
steeper increase with 6, related to the much higher value of V ■ Btest, and particularly so for Model 2. 

The location of the minimum of each of the orange curves is at dmin = -Zsmix/^^div, therefore its 
location depends on the average orientation and amplitude of the divergence field Bdiv with respect 
to the solenoidal field Bs. The orientation and amplitude of Bdiv also determines the height of the 
minimum (since the energy of the test field is fixed). With both divergence models, there are no 
general rules; i.e., the energy can increase or decrease with 6, and the location of the minimum 
depends on the case. 



7.3. Comparison witii tine potentiai fieid energy 

The physically meaningful quantity is represented by the energy normalized to the energy of the 
corresponding potential field, represented in Figure E]by black lines. For different 6 values, the 
normal component of the field Bdiv at the boundary changes according to Eq. (fTTt . hence also the 
energy of corresponding potential field depends — quadratically — on 6. Due to the additional 6- 
dependence, the shape of the black lines is not always parabolic in the six cases, and the actual 
profiles depend on the details of the spatial distribution of divergence in the test field. 

To show that, we first notice that the two divergence models behave very differently, except 
for Bdd where the range in 6 is too narrow to show significant differences. For instance, E/Ep of 
Model 1 (continuous black lines) is an increasing function of 6 in the range (-15, 15) in the Btd 
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case. Model 2, on the other hand, has a paraboUc energy profile with minimum atd ~ -4. For both 
models, the energy variation is relatively large (1.8 and above 2 for Models 1 and 2, respectively), 
whereas the variation in the same range of 5 is smaller for the Bdd and Bmhd cases. 

The extrapolated cases yield not only much larger variations (note again the difference in scales 
between the top and bottom rows of Figure |2]i, but also a stronger dependence on 6. In particular. 
Model 2 yields a relative energy that sharply increases with 6, for instance, to one order of mag- 
nitude increase for 6 going from the value to 1 in the Bex2pp case. A saturation at high values 
of 6 is clearly visible in the dashed black line (Model 2) of the Bexi case, and is hinted at in the 
Bex2pp case. Such saturation is actually present in all three extrapolated cases, yielding values that 
are higher than those shown in the corresponding plots. The saturation happens when the quadratic 
dependence on 6 of the energy of potential field compensates the quadratic term d'Eny. 

On the other hand. Model 1 shows a more complex dependence on 6, which is shown in mag- 
nified scale by the black lines in Figure |3] Counterintuitively, the largest variation in the relative 
energy E/Ep as a function of 6 is found for Bex2pp, i-e., for the extrapolation case, which satis- 
fies Thomson's theorem better, see Table [T] The continuous black lines in Figure |3}J,f show the 
presence of one maximum and one minimum in the considered range of values of 6 (for the Bex2pp 
case, these lie outside the considered range), implying that, at high values, the potential field energy 
grows faster than the total energy. The location of the extrema is different in the three Bexi , Bex2pp, 
and Bex2 cases, and in none of the cases are the extrema found for the solenoidal (6 - 0) or the test 
(6-1) configurations. In general, the maximum and minimum energy configurations depend on 
the spatial distribution of the divergence of the test field, through iSs.div 

7.4. Unphysical cases. 

The black dash-dotted line at E/Ep = 1 in Figure [3] is the value below which unphysical fields 
are obtained. We find that only Model 1 can produce unphysical solutions, and only for specific 
range of 6 values in the Bdd, Btd, Bex2 cases. The latter case is known from the value of Ep/E in 
Table [T] and is considered to be an extreme case because of the large divergence that it involves. 
However, the possibility of also creating unphysical solutions in the far more solenoidal field Bdd 
and Btd (for values of |<5| > 5) is unexpected. It confirms that not just the value of the divergence 
is important, but also its detailed spatial distribution with respect to the solenoidal component, as 
evident from /Ss.div It is the alignment between Bjiv and Bj, and not just the magnitude of Bdiv, 
which determines how strongly the energy depends on 6. Moreover, while E > Ep is always 
satisfied for Bexi and Bex2pp, the minimum value of E is close to Ep (see Figure |3}i,e), showing 
that unphysical fields may be found relatively easily in NLFFF extrapolations. 

From Table|2]and the related discussion of Sect. l5.3l we showed that the main source of violation 
of Thomson's theorem is the term fip.s/j.ns in ^mix- The dependence on 5 of this term, normalized 
to the energy of the test field, is shown by the red curves in Figure[3]for both models of divergence 
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Fig. 3. Magnetic energy, Eq. l ll6t . normalized to tiie energy of the corresponding potential field (black lines, 
zoom of Figure|2}, and energy associated to tiie nonsolenoidal part of Bj, normalized to the energy of Btcst (red 
lines, equal to Ej„JE in the notation of Eq. ((Sj), as a function of the amplitude 6 of the nonsolenoidal term, 
Eq. ( 114b . Two models for the divergence are shown according to Eq. l ll5t . Each panel shows the results of a 
test field: (a) the potential field of a double dipole, Bod\ (b) the TD model, Bxd; (c) the MHD model, Bmhd; 
(d) the NLFFF model of the nonpreprocessed magnetogram of AR 1 1 158, Bexi ; (e) the NLFFF model of the 
preprocessed magnetogram of AR 1 1024, Bex2pp; (f) the NLFFF model of the nonpreprocessed magnetogram 
of AR 11024, Bex2- The black dash-dotted line at E/Ep = 1 marks the value below which the solution is 
unphysical. A large change of scale of both axes is present between the top and the bottom rows. 



(Eq. (fTSll). The contribution to the total energy is negligible in the Bdd and Bmhd cases, and can 
be a few percent for large 6 in the Btd case. 

In the extrapolated cases, the dependence of iSp.s/j.ns on 6 is linear for Model 1 and parabolic 
for Model 2. In Model 1, the steepness of the linear curve increases, going to Bex2pp to Bexi and 
Bex2, as expected (see Table [T] and related text). The amplitude of the error is two to three orders 
of magnitude larger than in the Bdd, Bjd, and Bmhd cases. In the Bex2pp case the error is smaller, 
but it is stiU about a factor 20 larger than in Btd for 6 = 5. 

If we consider the black curves in Figure |3] for 5 = 0, we can identify the energy of the 
solenoidal field as a natural reference value for the free energy. Starting from this reference value, 
for increasing \6\, the linear contribution of /ip^s/j.ns, together with the quadratic change in the po- 
tential field energy, creates the maximum and minimum values of E/Ep. If the linear contribution 
is large enough, the minimum lies below the threshold E/Ep = 1, and there is a range of values 
where the solution is unphysical. For even higher values of \6\, the quadratic dependence of the 
potential field energy dominates E. From this point onward, £p,s/j,ns is not the main source of error 
in Eq. ^. 

More generally, a parametric study like the one in Figs. |2]and|3]can be used to identify what 
is the level of divergence (i.e., the level of fimixOr £p,s/j,ns) that can be tolerated and which is the 
threshold above which the solution becomes entirely unphysical [i.e., with E/Ep < 1). 
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Fig. 4. Power spectra of the two-dimensional fields (continuous line) and B- (dashed line), summed over 

all wave numbers ky, for the three cases Btd (orange), Bex2pp (red), and BekI (black). Left: at the bottom 
boundary. Right: at the tenth pixel in height. Spectra are normalized to their maximum value, and the spatial 
resolution is taken to be 1 in both directions (i.e., the wave number fc, has the dimension of pixel"' and is 
normalized to the total number of modes). 

In conclusion, the parametric study shows that the energy may be severely influenced by the 
solenoidal property of the field. The effect depends not only on the amplitude of the nonsolenoidal 
component, but also on the specific average orientation of the nonsolenoidal component with re- 
spect to the solenoidal one (directly affecting Zss div in Eq- (fTSIl). As a result, a single-number 
divergence metric, such as ( ), is insuflicient to deduce what errors should be expected in the 
energy. A more proper indication is found by the numerical verification of Thomson's theorem 
(Sect. |5]i and by a parametric study as presented in this section. 

8. Source of divergence in NLFFF extrapolations 

We now investigate in more detail some of the test fields discussed in Sect. |5] with emphasis on the 
reason for the large divergence that leads to violating Thomson's theorem. The main source of error 
comes, in almost all the cases, from the mixed term fimix, and is associated with the nonsolenoidal 
component of the current-carrying part of the field. Also, there are markedly larger errors in the 
extrapolated test fields, Bexi, Bex2pp, Bex2, than in Bdd, Btd, and Bmhd- Finally, the preprocessing 
of the vector magnetogram before extrapolation yields more solenoidal fields, whereas a simple 
averaging does not seem to be enough for removing errors, and yields a more severe violation of 
Thomson's theorem (Eq. 

8.1. Analysis of small scales 

One main difference among the Btest cases in the upper half of Table [1] is the length scale of the 
magnetic field: While the first three cases are smooth fields with a magnetic field variation spanning 
several times the spatial resolution, the extrapolated cases have large variations on the pixel scale, 
especially at the bottom boundary, i.e., on the vector magnetogram that is used as a boundary 
condition for extrapolations. This is true to a different degree for the three cases: For Bexi the 
vector magnetogram was interpolated (with a flux-conserving average) at a resolution of about one 
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third that of Bex2 and Bex2pp- Such an interpolation smooths part of the small scale away, yielding 
results that are closer to the Bex2pp case rather than to the Bex2 one. Bex2pp is not interpolated, but 
it is preprocessed, an operation that includes an explicit smoothing of smaller scales, especially on 
the transverse components. Finally, Bex2 has neither interpolation nor preprocessing, and it retains 
all the small scales that are present at the full resolution of the Hinode/SOT vector magnetograms. 

As an example, Fig.|4]shows the power spectrum of the x- and z-components of the fields Btd, 
Bex2pp, and Bex2, at two different heights as a function of the normalized wave number k^ . The 
lefthand panel of the figure shows that, at the bottom boundary, Bjd has power spectra that decrease 
rapidly with k^, in both components. In contrast, the power spectra of Bex2pp and Bex2 have higher 
values on all scales, which are particularly strong in the vertical component. 

Ten pixels above the bottom boundary (right panel in Fig. HJ, the Btd power spectrum is es- 
sentially the same as at z = because both planes cut through the flux rope, so a similar magnetic 
structure is present. In contrast, Bex2pp on the upper plane has a much more peaked spectrum, 
except for the distribution tail on the smallest scales which is basically as strong as at the bottom 
boundary. Such a component on the shortest scales comes from the force-free condition that is en- 
forced by the extrapolation code, which propagates into the volume the small scales that are present 
at the bottom boundary. 

We now consider the difference between preprocessed case Bex2pp and the non-preprocessed 
one Bex2- The difference in { \fi\ ) between the two is about a factor 2, and it is large in the other 
energy metrics in Table [T] The comparison between the normalized spectra of Bex2pp and Bex2 in 
Fig. m shows that there are comparable (relative) energies on small scales in both cases. Actually, 
by locally changing the magnetic field at the bottom boundary to enforce force-free compatibility, 
preprocessing increases the small scales. The smoothing term that is present in preprocessing only 
has a limiting effect on such an increase. Therefore, the two cases Bex2pp and Bex2 do not differ 
strongly as far as the presence of small scales is concerned, while Thomson's theorem is much 
better satisfied for Bex2pp than for Bex2 (see Sect. 15. lb . 

The cleaned test fields Btest.s are numerically solenoidal, and there is no violation of Thomson's 
theorem. However, in these cases, too, small scales are increased (not shown), since the divergence 
cleaner introduces extra electric currents that are related to derivatives of the divergence of the 
original field, see Eq. (IB. 5b . This is an additional confirmation that the presence of small scales as 
such is not directly at the origin of the violation of Thomson's theorem. 



8.2. Role of small scales and preprocessing 



Valori et al 



( 120101) show that the NLFFF extrapolation of the Bjd vector magnetogram yields a 
very accurate reconstruction of the whole test field, which is also solenoidal to a very high degree. 
On the other hand, there is a large difference in the scale distribution between smooth fields like 
the Btd and the extrapolated fields. 
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The presence of small scales inside the volume, which are induced by the small scales at the 
boundary, may not be correctly approximated by the discretization employed in extrapolation code, 
yielding local violation of the solenoidal constrain. However, when the extrapolation from a pre- 
processed magnetogram is considered, the extent of the violation of Thomson's theorem is greatly 
reduced, even though small scales are actually increased. By partially enforcing force-free compat- 
ibility on the bottom boundary, the preprocessing provides the extrapolation code with a boundary 
condition that is more compatible with the force-free equations. Since extrapolation codes attempt 
to construct a solution of the force-free equations that is simultaneously force- and divergence- 
free, the more compatible the boundary, the more consistent (i.e., force- and divergence-free) the 
obtained solution. Conversely, when the boundary condition is incompatible with the force-free 
equation, the reduction of the Lorentz forces is at the expense of the solenoidal condition. In such 
cases, the divergence of the solution is higher, and Thomson's theorem is more severely violated. 
We thus conclude that the incompatibility of the boundary condition with the force-free condition 
is at the origin of the difference in the errors Zsj „, and ^mix between Bex2pp and Bex2- 

We notice that preprocessing is a parametric method that can produce progressively more force- 
free-compatible vector magnetograms for higher values of the employed parameters, at the price 
of larger modifications of observed values. The energy values and their relative errors therefore 
vary continuously as a function of the preprocessing pa rameters, quite independently of the partic - 



ular extrapolation method that is employed (see, e.g.. 



Schriiver et al 



2008 



Metcalf et al. 



2008h . 



No unequivocal method is available in order for determining the best pa rameters to use (see, e.g.. 



Wiegelmann et al 



2006 



Fuhrmann et al 



2011 



Wiegeknann et al 



20121), which leaves energy es- 



timations subjected to uncomfortable arbitrariness. 



9. Conclusions 

Thomson's theorem states that the energy of a magnetic field is given by the sum of the energy of the 
current-carrying part of the field plus the energy of the potential field that has the same distribution 
of the normal component on the boundary of the considered volume. The field must be perfectly 
solenoidal for the theorem to be valid. Such a condition is often only approximately satisfied in 
numerical simulations, such as in MHD simulations and NLFFF extrapolations. However, it is 
a non-trivial task to identify a quantitative estimation of solenoidal errors that can be applied to 
different discretizations of magnetic fields, essentially due to the non-local consequences that such 
errors produce. Our goal has been to develop physically meaningful metrics and practical methods 
that can be used to judge whether the solenoidal property is fulfilled with sufficient accuracy. 

To this aim, we introduced a decomposition of the energy of a discretized field into solenoidal 
and nonsolenoidal contributions that allowed an unambiguous and numerically well-defined esti- 
mation of the effect of the divergence in terms of associated energies. Moreover, we introduced a 
method of parametrizing the divergence that allows for an exploration of the nonsolenoidal effects. 
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In this way, the numerical verification of Thomson's theorem offers an operational and quanti- 
tative way of checking the rehabihty of energy estimations in numerical computations. Since the 
violation of Thomson's theorem is solely determined by the presence of magnetic charges, it is at 
the same time a quantitative estimation of the importance of solenoidal errors. 

We applied our method to six different test cases, covering a representative sample of numerical 
realizations. Of the six test cases considered here, two of them (the dipolar field Bdd and a snapshot 
of an MHD simulation of null-point reconnection Bmhd) presented negligible violations, and one 
(a force-free current ring Btd) offered only a moderate one that, however, has finite effects on the 
energy. In the case of an NLFFF extrapolation of a preprocessed vector magnetogram (Bex2pp), 
the sum of the potential energy Ep and free energy £j_s is very close to the total energy E, and 
one could draw the conclusion that almost no violation of Thomson's theorem occurs. However, 
by separating all contributions in Eq. (|7), our analysis showed compensating energy contributions 
(iij^ns and iiniix) that are close to Ej ^- If the most consei-vative view is adopted by considering errors 
in absolute values, then the opposite conclusion must be drawn: The violation is large enough to 
compromise the estimation of the free energy, since both Ej ^^ and iSmix are on the order of the free 
energy value Zsj,,. The last two cases we studied, also NLFFF extrapolations but of nonpreprocessed 
magnetograms (Bexi and Bex2), represent cases with very large errors. 

The energy of the potential field , is the reference value for the free energy. In our appli- 
cations, the inaccuracy in its determination, Ep^s, is practically never significant. The current- 
carrying part of the field is responsible for the largest errors instead. 

The parametric study shows that the amplitude of the nonsolenoidal component is not the only 
factor that generates errors in the energy. The average orientation of the nonsolenoidal component 
with respect to the solenoidal one (affecting directly fis.div in Eq. ( fTSI l) plays an even more important 
role. Indeed, even using a relatively solenoidal discretized magnetic field (like Btd), it is possible to 
create configurations where the energy of the field is lower than that of the corresponding potential 
field. Such unphysical solutions have also been found in some cases of NLFFF extrapolations. 

More generally, in NLFFF extrapolations the energy of the reconstructed field was found to 
vary according to the extent of the modification that was enforced on the vector magnetogram that 
is used as boundary condition (by preprocessing, i.e., by smoothing and/or by enforcing force- 
freeness compatibility). Our study shows quantitatively the effect of such practices on the energy, 
and makes it clear that the origin of the variability (and errors) in energy estimations based on 
NLFFF extrapolations is the presence of a large divergence, which is eventually caused by the lack 
of compatibility between the equations solved (solenoidal force-free field) and the photospheric 
boundary conditions, rather than by noise or the small scales present in the vector magnetogram. 

Finally, the parametric study is based on a numerically solenoidal field that is obtained from a 
given, nonsolenoidal one. We introduced a method for the complete removal of the nonsolenoidal 
component of a discretized field. At the price of changing boundary values and the current density, 
this method provides a field that is solenoidal to numerical precision. When the solenoidal versions 
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of the test fields are considered, the Thomson theorem is found to be fulfilled with more than 99% 
accuracy. 

We concluded that testing Thomson's theorem in numerical realizations of magnetic fields is a 
powerful method quantifying the amount of nonsolenoidal contributions to a numerical magnetic 
field. In particular, it allows assessing the reUabiUty of free magnetic energy estimations, a crucial 
quantity in phenomena such as flares and coronal mass ejections. To this purpose, we proposed 
a set of analytical and numerical tools that allowed us to fully test the reliability of numerical 
magnetic fields. Such a set includes a method for removing the divergence from a given discretized 
field, to numerical precision. The effect of larger and larger divergence contributions is studied 
by parametrically adding a known divergence to the numerically solenoidal field. In this way, it 
is possible to monitor the effect of the nonsolenoidal part of the magnetic field and to quantify its 
effect in terms of magnetic energy. Our method can be applied to any discretization of magnetic 
fields, e.g., in MHD simulations and in NLFFF extrapolations, to constrain quantitatively errors 
due to violation of the solenoidal property. 

Acknowledgements. The authors are pleased to thank Guillaume Aulanier for fruitful discussions, Bemhard Kliem and Ti- 
bor Torok for making the numerical solution of the TD equilibrium available, and the referee for helpful comments that 
improved the clarity of the paper. GV is indebted to the NLFFF Consortium for stimulating discussions and collaborations. 
The research leading to these results has received funding from the European Commission's Seventh Framework Program 
(FP7/2007-2013) under the grant agreement eHeroes (project n° 284461, www.eheroes.eu). SM gratefully acknowledges 
support from the NASA Postdoctoral Program, administrated by Oak Ridge Associated University through a contract with 
NASA, during her stay at NASA Goddard Space Flight Center Calculations were done on the quadri-core bi-Xeon comput- 
ers of the Cluster of the Division Infonnatique de I'Obsei'vatoire de Paris (DIO). SDO data are courtesy of NASA/SDO and 
the HMI science team. Hinode is a Japanese mission developed and launched by ISAS/JAXA, collaborating with NAOJ as a 
domestic partner, NASA and STFC (UK) as international partners. Scientific operation of the Hinode mission is conducted 
by the Hinode science team organized at ISAS/JAXA. This team mainly consists of scientists from institutes in the partner 
countries. Support for the post-launch operation is provided by JAXA and NAOJ (Japan), STFC (UK), NASA (USA), ESA, 
and NSC (Norway). 

References 

Amari, T., Luciani, J. R, Aly, J. J., Mikic, Z., & Linker, J. 2003, ApJ, 585, 1073 

Asai, A., Ishii, T. T, Kurokawa, H., Yokoyama, T, & Shimojo, M. 2003, ApJ, 586, 624 

Aulanier, G., Demoulin, P, Schrijver, C, et al. 2013, A&A, 549, A66 

Aulanier, G., Janvier, M., & Schmieder, B. 2012, A&A, 543, Al 10 

Aulanier, G., Torok, T., Demoulin, P, & DeLuca, E. E. 2010, ApJ, 708, 314 

Berger, M. A. 1985, ApJS, 59, 433 

Brackbill, J. U. & Barnes, D. C. 1980, Journal of Computational Physics, 35, 426 

Chen, P. F. 20 11, Living Reviews in Solar Physics, 8, 1 

Emslie, A. G., Dennis, B. R., Shih, A. Y., et al. 2012, ApJ, 759, 71 

Fan, Y. 2010, ApJ, 719, 728 

Forbes, T. G. 2000, J. Geophys. Res., 105, 23153 

Forbes, T. G., Linker, J. A., Chen, J., et al. 2006, Space Sci. Rev, 123, 251 

Fuhrmann, M., Seehafer, N., & Valori, G. 2007, A&A, 476, 349 

Fuhrmann, M., Seehafer, N., Valori, G., & Wiegelmann, T. 2011, A&A, 526, A70 

Golub, L. & Pasachoff, J. M. 2009, The Solar Corona 

Karpen, J. T, Antiochos, S. K., & DeVore, C. R. 2012, ApJ, 760, 81 

Lawrence, C. E. 1998, Partial Differential Equations (American Mathematical Society) 

Liu, R., Liu, C, Wang, S., Deng, N., Wang, H. 2010, ApJ, 725, L84 

Low, B. C. 2001, J. Geophys. Res., 106, 25141 

Article number, page 26 of 129 1 



G. Valori et al.: Accuracy of magnetic energy computations 

Masson, S., Aulanier, G., Pariat, E., & Klein, K.-L. 2012, Sol. Phys., 276, 199 
Metcalf, T. R., Derosa, M. L., Schrijver, C. J., et al. 2008, Sol. Phys., 247, 269 
Priest, E. R. 2003, Solar magnetohydrodynamics, ed. B. N. Dwivedi, 217-237 
Savage, S. L., McKenzie, D. E., Reeves, K. K. 2012, ApJ, 747, L40 
Schrijver, C. J., DeRosa, M. L., Metcalf, T., et al. 2008, ApJ, 675, 1637 
Schrijver, C. J. & Zwaan, C. 2008, Solar and Stellar Magnetic Activity 
Shibata, K. Sl Magara, T. 201 1, Living Reviews in Solar Physics, 8, 6 
Torok, T. & Kliem, B. 2003, A&A, 406, 1043 
Tajima, T. & Shibata, K. 2002, Plasma astrophysics 
Taylor, J. B. 1986, Reviews of Modem Physics, 58, 741 
Titov, V. S. & Demoulin, P 1999, A&A, 351, 707 
Valori, G., Demoulin, P, & Paiiat, E. 2012, Sol. Phys., 278, 347 
Valori, G., Green, L. M., Demoulin, P, et al. 2011, Sol. Phys., 374 
Valori, G., Kliem, B., Torok, T., & Titov, V. S. 2010, A&A, 519, A44 
Wheatland, M. S., Sturrock, P A., & Roumeliotis, G. 2000, ApJ, 540, 1150 
Wiegelmann, T., Inhester, B., & Sakurai, T. 2006, Sol. Phys., 233, 215 
Wiegelmann, T., Thalmann, J. K., Inhester, B., et al. 2012, Sol. Phys., 281, 37 

List of Objects 

Appendix A: Details of the numerical implementation 

In the applications presented in this paper we have considered uniform Cartesian grids of resolution 
A in all directions, discretizing a rectangular volume 'V (see Sect.|4]for the actual values of A and 'V 
in each case). We compute derivatives using the standard second-order, central-difference operator, 
and we employ the relevant one-sided {i.e., forward or backward), second-order differences at the 
boundaries of "V. The only exception is the computation of the divergence of Btest, since all test 
fields are known in a volume that is larger than the selected (on lateral and top boundaries). In 
this case, V • Btest is computed using the central differences also at the location of the lateral and 
top boundaries of "V. 

In the computation of volume integrals, the cell volume A^^ is assigned to each internal node 
of the grid, whereas the cell volume is reduce to half, one fourth, and one eighth for nodes on the 
lateral surfaces, edges, and corners of "V, respectively. Similarly, in the computation of surface in- 
tegrals, the cell surface A^ is assigned to each node inside each side of "V, whereas the cell surface 
is reduced to half and one fourth on edges and comers of each side, respectively. Despite the accu- 
rate computation of integrals, the divergence theorem, Eq. (|4]i, is not insured to hold numerically, a 
property that requires special techniques, like finite-volume discretizations, to be fulfilled. 

Appendix B: Divergence cleaner 

To construct a numerically solenoidal field [Bg] from a field [B] let us define 



Bs = V X A, 



(B.l) 
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where A is the vector potential computed from B in the volum e = [^^1,^1:2] x [yi,3'2] x Ui,22]- 



Valori et al 



(120121) using the gauge z ■ A = 0, yielding 



The vector potential A can be derived as in 
the expression 

A = b + zx^'fidz', (B.2) 
where b = {A^{x, y,z = Z2), Ay{x, y,z = Z2), 0) is any solution of 

= d,by - dyb, - B„,,,(x, y, z = Z2) . (B.3) 
A direct substitution of Eq. (IB.2b into Eq. (IB. lb shows that 

s V X A = B + z J" '(V ■ B) dz' , (B.4) 



with the property that V ■ Bj = 0. In other words, Eq. (IB .4b naturally separates B into a solenoidal 
part Bj and a nonsolenoidal one, thus defining a divergence cleaner for B. The z-component of B is 
changed throughout the volume except on the top boundary, whereas the x- and y-components are 
unchanged. The amplitude of the modification to B at a given height z is given by the cumulative 
effect of "magnetic charges" above that altitude. 

Since only the z-component of the field is changed, the divergence cleaner changes the x- and 
y-components of the current, but not the z-component, 

Js = J + {dy, -5,, 0) J"(V • B) dz' , (B.5) 



therefore the cleaner changes the injected magnetic flux but not the injected electric current through 
the bottom layer. On the other hand, since most of the test fields considered in this article have 
the highest values of divergence close to the bottom boundary, only the lower part of the field is 
changed significantly by the cleaner. 

Computation Bj requires numerical computation of an integral of the type G(z) - J_ ' f{t) df, 
as in Eq. ( IB.4I ) for / = V • B. To achieve numerical accuracy in the solenoidal property of B,, 
G(z) must satisfy dfi{z) - -f{z) numerically, i.e., must satisfy the numerical formulation of the 
fundamental theorem of integral calculus in the employed discretization. For the second-order 
central differences that are used in the analysis, this can be obtained by the recurrence formulae 

G(«z-1) = 0, 

G{k) = G(/t + 2)-H2A/(yt-i- 1), < /t < «z - 3, (B.6) 

where G(z) = G(zi +A:A) = G{k) with A: = 0, 1 , 2, ■ ■ ■ , (n^- 1), and A is the uniform spatial resolution 
in z. 

The constraint d,G(z) = -f{z) in the second-order, central-difference discretization does not 
fix the value of G(nz - 2). To do that, we require that the divergence of Eq. (IB .4b also vanishes at 
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the bottom boundary, i.e., (V ■ Bs)|z=zi - 0. Here the second-order divergence operator is computed 
by using a second-order, forward derivative in the z-direction, i.e., defining the operator Vqs = 
V,,y + zdf, where V,,,, = xd, + fdy and (5f/)(0) - (-3/(0) + 4/(1) - /(2))/2A. By using 
the recurrence formula Eq. (IB. 6) , the condition on the bottom boundary is transformed into the 
condition for G(nz - 2), yielding 

m+ ^[/os(0)-3/(i)] , 

where /os - Vqs ■ B. Such a numerical trick is only possible if the volume is discretized by an even 
number of points in the z-direction, therefore the analysis volumes employed in the article were 
chosen to satisfy such a requirement. 



G{n, - 2) = A 



2 - Z 



Appendix C: Measures of V • 6 



The total divergence of a field B can be conveniently expressed by a single number using the 
average ( |/i) over the grid nodes of the fractional flux 



jdv (V-B),. 



(C.I) 



20001) . Taking a 



through the surface dv of a small volume v including the node ; jWheatland et al. 
cubic voxel of side equal to A as the small volume v centered on each node, the divergence in the 
discretized volume of uniform and homogeneous resolution A is then given by 



<l/l) 



A 

6N 



V • B, | 



4^ IB/ 



(C.2) 



where i runs over all nodes in 'V. This metric depends on the considered volume, so that values 
are strictly comparable only if computed on equal volumes. 
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